Univariate visualization

library(reshape2)
# Let's take the average of the three results
lakes_complete <- (imp1 + imp2 + imp3)/3.0
lakes_complete$URBAN <- imp1$URBAN
lakes_complete$LAKE_ORIGIN <- imp1$LAKE_ORIGIN

final_lakes <- cbind(UID=select_lakes2$UID, lakes_complete)
#write.csv(final_lakes, "final.raw.12.csv")

final_summ <- summaries(lakes_complete[,1:43])

final_summ %>% kable()%>% kable_styling(bootstrap_options = "striped", full_width = F, font_size = 12)

lakes_long <- melt(final_lakes[,1:44], id.var = "UID")
head(lakes_long)

##       UID variable value
## 1 1000001 ALUMINUM 0.007
## 2 1000010 ALUMINUM 0.020
## 3 1000011 ALUMINUM 0.012
## 4 1000013 ALUMINUM 0.031
## 5 1000014 ALUMINUM 0.000
## 6 1000015 ALUMINUM 2.480

ggplot(lakes_long, aes(x=value))+
  geom_histogram(aes(y=..density.., fill="chartreuse4", color="chartreuse4"), bins = 20)+
  geom_density(fill="deepskyblue", alpha=0.2)+
  facet_wrap(~variable, ncol = 5, scales = "free")+
  theme_bw()+guides(fill=F, color=F)

From the above plot we can see that most of the histograms are highly right-skewed with most of the data having very low value and few high values. This is normal in such kind of environmental systems indicating most of the lakes have low baseline concentration and few of the lakes have high values that could be an indication of human influence or the characteristic of that specific lake. To make the data more or less similar to normal distribution we can log transform the skewed features.

In addition, some of the features have value of zero this could be that it is actually zero or could be as a result of detection limit. The EPA reports the detection limit for each feature, the detection limit for each feature varies. In order to transform we want the features to be strickly positive, and for this we will add the smallest detection limit of all the feature. As noted in the summary section the elevation and ANC features have some negative values, to deal with this I will use cube root transformation on them before using log transformation.

library(dplyr)

# Add small purturbation to the data to avoid the zeros in the data
small_c <- 0.002
final_lakes.c <- final_lakes[,2:44] + (final_lakes[,2:44] < 0.002)*small_c

# Check if this has a significant effect on the data
new_summ <- summaries(final_lakes.c)
sum_new <-  new_summ - summaries(final_lakes[,2:44])
head(sum_new,10)

##             min max sd  mean NAs
## ALUMINUM  0.002   0  0 0.001   0
## AMMONIA_N 0.002   0  0 0.000   0
## ANC       0.002   0  0 0.000   0
## CALCIUM   0.000   0  0 0.000   0
## CHLORIDE  0.000   0  0 0.000   0
## COLOR     0.002   0  0 0.000   0
## COND      0.000   0  0 0.000   0
## DOC       0.000   0  0 0.000   0
## MAGNESIUM 0.000   0  0 0.000   0
## NITRATE_N 0.002   0  0 0.001   0

After inspecting the difference between the pruturbed and the new data, there is a slight change which is in acceptable range. Now we can start transforming the data either using log or cube root, whichever works on the specific feature. Looking at the data two features, Elevation and ANC (acid neutralizing capacity), have some negative values, and we can apply log transformation directly. Either we have to shift the data to be stricktly positive or we can try cube root transformation. First, I will try to use cube root transformation on these features, and if the result is not satisfactory, we can try to shift and log transform.

# Copy the data frame
transform_lake.df <- final_lakes.c
# Cube root transform first
cube.r <- function(x){sign(x)*abs(x)^(1/3)}
transform_lake.df$ANC <- cube.r(transform_lake.df$ANC)
transform_lake.df$ELEVATION <- cube.r(transform_lake.df$ELEVATION)

Check their histogram before and after transformation.

par(mfrow=c(2,2))
hist(final_lakes.c$ELEVATION, breaks=50, main="Original", col="#18BF5C", freq=F, xlab="Elevation")
lines(density(final_lakes.c$ELEVATION), col="darkred", lwd=3); box()
hist(transform_lake.df$ELEVATION, main="Transformed", breaks=50, col="#18BF5C", freq=F, xlab="Elevation")
lines(density(transform_lake.df$ELEVATION), col="darkred", lwd=3); box()

hist(final_lakes.c$ANC, breaks=50, main="Original", col="#18BF5C", xlab="ANC"); box()
lines(density(final_lakes.c$ANC), col="darkred", lwd=3); box()
hist(transform_lake.df$ANC, main="Transformed", breaks=50, col="#18BF5C", freq=F, xlab="ANC")
lines(density(transform_lake.df$ANC), col="darkred", lwd=3); box()

As it can be seen in the above histograms the elevation already got transformed effectively to gaussian distribution, while the ANC might still needs additional transformation, but I will leave it as it is.

# Variables that don't need transformation
ok_vars <- c("ANC","PH", "OXYGEN", "ELEVATION","INDEX_LAT_DD", "INDEX_LON_DD", "MMI_BENT_NLA12", "MMI_ZOOP_NLA6", "RDis_IX", "URBAN", "LAKE_ORIGIN")

# Variables that need transformation
log_vars <- c("ALUMINUM","AMMONIA_N", "CALCIUM", "CHLORIDE", "COLOR", "COND", "DOC", "MAGNESIUM", "NITRATE_N","NITRITE_N","NTL", "POTASSIUM", "PTL","SILICA", "SODIUM", "SULFATE","TURB", "CHLX", "CHLL", "TEMPERATURE", "DEPTH", "CONDUCTIVITY", "SECCHI", "AREA_HA", "ABUNDANCE", "PERIM_KM", "AMITOTAL", "LitCvrQc3OE", "LitRipCvrQc3OE", "RVegQc3OE", "MICX","ATRAZINE","METHYLMERCURY", "TOTALHG")

#for (i in 1:length(log_vars)) {
#  transform_lake.df[,log_vars[i]] <- log(transform_lake.df[,log_vars[i]])
#}

transform_lake.df[, log_vars] <- log(dplyr::select(transform_lake.df, log_vars))

transform_lake.df$UID <- final_lakes$UID
#trans_summ <- summaries(transform_lake.df)

ts_lake_long <- melt(transform_lake.df, id.var = "UID")

ggplot(ts_lake_long, aes(x=value))+
  geom_histogram(aes(y=..density.., fill="#81D2FF", color="#307D00"), bins = 10)+
  geom_density(fill="#FFB05A", alpha=0.2)+
  facet_wrap(~variable, ncol = 5, scales = "free")+
  guides(fill=F, color=F) + theme_bw()

As it can be seen, most of the features have a unimodal normal distribution except few. For now, I am going to leave it as it is and proceed with the next section.

Bivariate plot

library(corrplot)

## corrplot 0.84 loaded

# Compute correlation table
lakes_cor <- cor(transform_lake.df[,-44])
#dev.new(width=10, height=10, unit='in')
clrs <- colorRampPalette(c("#B26205","#FFB05A", "white","#81D2FF", "#0980B2"))
corrplot(lakes_cor, type = "lower", tl.pos = "ld", mar = c(1,1,1,1), 
         tl.col="black", addCoef.col = "black", tl.cex = 0.9, 
         diag = F, method = "square",col=clrs(100), number.cex = 0.55)

Next, let us filter highly correlated features and visualize them.

# Filter those with high correlation
threshold <- 0.85
high_cor <- lakes_cor
diag(high_cor) <- 0
ok <- apply(abs(high_cor) >= threshold, 1, any)
high_cor <-lakes_cor[ok,ok]

corrplot(high_cor, type = "lower",tl.cex = 1.2, tl.pos = "ld", mar = c(1,1,1,1), 
         addCoef.col = "black", method = "square", number.cex = 1, 
         tl.col="black", col=clrs(100), diag = F)

Most of the features having high correlations are expected fron what we know about water quality parameters. For example, CHLL is littoral chlorophyl-a concentration and CHLX is chlorophyl-a concentration in the lake, and this two are expected to have high correlation. In addition, if we look at COND (condactivity of lake’s water) which is due to the ions present in water such as calcium. chloride, sodium, etc., we expect high correlation with them. A more detailed plot is also shown below.

library(psych)
high_cor_names <- colnames(high_cor) 
high_cor_lakes <- subset(transform_lake.df,select = high_cor_names)

pairs.panels(high_cor_lakes, 
             main="Pairs plot for highly correlated features",
             method = "pearson", # correlation method
             hist.col = "#0980B2",
             density = TRUE,  # show density plots
             ellipses = TRUE # show correlation ellipses
             )

Map plot using plotly

We can also see the exact locations of the the sampled lakes as shown below. We can zoom and hover on the points to see the total nitrogen, total phsophorus and area of the lakes.

library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

#final_lakes.ts <- transform_lake.df
# geo styling
g <- list(
  scope = 'usa',
  projection = list(type = 'albers usa'),
  showland = TRUE,
  showlakes = TRUE,
  landcolor = toRGB("gray95"),
  subunitcolor = toRGB("gray85"),
  countrycolor = toRGB("gray85"),
  countrywidth = 0.5,
  subunitwidth = 0.5
)


p <- plot_geo(final_lakes, lat = ~INDEX_LAT_DD, lon = ~INDEX_LON_DD, color= ~CHLX) %>%
  add_markers(
    text = ~paste(paste("Total nitrogen:",round(NTL,2)), paste("Total phosphorus:",round(PTL,2)), paste("Area (HA):", round(AREA_HA,2)),sep = "<br />"),
    symbol = I("square"), size = ~AREA_HA, hoverinfo = "text"
  ) %>%
  colorbar(title = "Chlorophyl a ammount<br />In the lake") %>%
  layout(
    title = 'US lakes<br />Surveyed in 2012', geo = g
  )

## Warning: `line.width` does not currently support multiple values.

p

Evaluation of multiple linear regression on the test data

To know how strongly the linear model trained above performed, we will use the test data and evualate its performance. Later, we will use this to compare it with other models.

performance <- function(y.obs, y.pred){
   model.rmse <- RMSE(y.obs, y.pred)
   model.cor <- cor(y.obs, y.pred)
   model.mae <- mean(abs(y.obs - y.pred))
   model.r2 <- 1- sum((y.obs - y.pred)^2)/sum((y.obs - mean(y.obs))^2)
   return(rbind(model.rmse, model.cor, model.mae, model.r2))
 }

y.pred <- predict(chl.lm, lake_test )

# Models performance metrics
ml.perf <- performance(lake_test$Chlorophyl, y.pred)
ml.perf

##                 [,1]
## model.rmse 1.3691218
## model.cor  0.8915594
## model.mae  0.9881585
## model.r2   0.7944393

The model performance as expected decreased slightly on the test data, however, we still get fair performance with good correlation between the test and predicted values. The model gives a mean absolute error of 0.9881585, which corresponds to an error of \({e}^{chlorophll}=\) 2.686283 ug/L in chlorophyll concentration. Taking into consideration the range of chlorphlyll concentration, which is from 0.4406993 to 1166, this error can be considered acceptable.

Stepwise linear model selection

From the above multiple linear regression result we have observed that only 15 features had a significant linear relationship with the response variable. Hence, we can further narrow down the number of features and select only the important features for the linear model using the step function to improve it further. I will use the step function with both backward and forward selection to the feature selection.

# Define a base model - intercept only
base.mod <- lm(Chlorophyl ~ 1, data=lake_train)
# Define the full model - including all predictors
all.mod <- lm(Chlorophyl ~ ., data=lake_train)

lakes_step <- step(base.mod, scope = list(lower = base.mod, upper = all.mod), direction = 'both', k=2, trace = F)
summary(lakes_step)

## 
## Call:
## lm(formula = Chlorophyl ~ SECCHI + NTL + ABUNDANCE + TURB + ELEVATION + 
##     NITRATE_N + DOC + PTL + SILICA + MMI_BENT_NLA12 + COND + 
##     CHLORIDE + CALCIUM + ALUMINUM + AMITOTAL + RVegQc3OE + LAKE_ORIGIN + 
##     SODIUM, data = lake_train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.4683 -0.7674  0.0238  0.8267  5.7185 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         1.391495   0.752253   1.850 0.064659 .  
## SECCHI             -0.501315   0.091961  -5.451 6.38e-08 ***
## NTL                 1.256190   0.116888  10.747  < 2e-16 ***
## ABUNDANCE           0.568734   0.061356   9.269  < 2e-16 ***
## TURB                0.466093   0.062380   7.472 1.80e-13 ***
## ELEVATION          -0.083296   0.018339  -4.542 6.29e-06 ***
## NITRATE_N          -0.226141   0.032643  -6.928 7.92e-12 ***
## DOC                -0.324669   0.109876  -2.955 0.003206 ** 
## PTL                 0.338990   0.072154   4.698 3.02e-06 ***
## SILICA             -0.094126   0.040224  -2.340 0.019489 *  
## MMI_BENT_NLA12     -0.009974   0.003515  -2.838 0.004637 ** 
## COND               -0.594700   0.135480  -4.390 1.26e-05 ***
## CHLORIDE            0.241610   0.053448   4.520 6.96e-06 ***
## CALCIUM             0.163125   0.080454   2.028 0.042887 *  
## ALUMINUM           -0.130094   0.039055  -3.331 0.000899 ***
## AMITOTAL           -0.057223   0.028024  -2.042 0.041438 *  
## RVegQc3OE           0.081985   0.038706   2.118 0.034424 *  
## LAKE_ORIGINNATURAL -0.213220   0.099790  -2.137 0.032880 *  
## SODIUM             -0.136208   0.083607  -1.629 0.103617    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.364 on 943 degrees of freedom
## Multiple R-squared:  0.783,  Adjusted R-squared:  0.7789 
## F-statistic: 189.1 on 18 and 943 DF,  p-value: < 2.2e-16

lakes_step

## 
## Call:
## lm(formula = Chlorophyl ~ SECCHI + NTL + ABUNDANCE + TURB + ELEVATION + 
##     NITRATE_N + DOC + PTL + SILICA + MMI_BENT_NLA12 + COND + 
##     CHLORIDE + CALCIUM + ALUMINUM + AMITOTAL + RVegQc3OE + LAKE_ORIGIN + 
##     SODIUM, data = lake_train)
## 
## Coefficients:
##        (Intercept)              SECCHI                 NTL  
##           1.391495           -0.501315            1.256190  
##          ABUNDANCE                TURB           ELEVATION  
##           0.568734            0.466093           -0.083296  
##          NITRATE_N                 DOC                 PTL  
##          -0.226141           -0.324669            0.338990  
##             SILICA      MMI_BENT_NLA12                COND  
##          -0.094126           -0.009974           -0.594700  
##           CHLORIDE             CALCIUM            ALUMINUM  
##           0.241610            0.163125           -0.130094  
##           AMITOTAL           RVegQc3OE  LAKE_ORIGINNATURAL  
##          -0.057223            0.081985           -0.213220  
##             SODIUM  
##          -0.136208

We can see that step wise linear model selection has selected 18 features/variables as relevant to the linear model, which by itself is high dimension. Next, we can calculate the importance ranking of these features and visualize it.

#library(mlbench)

# estimate variable importance
predStepwise <- varImp(lakes_step, scale=FALSE)

#Change names
rownames(predStepwise)[11] <- "URBAN"
rownames(predStepwise)[14] <- "LAKE_ORIGIN"

par(mai=c(1.5,1,0.2,0.1))
bar.df <- predStepwise[order(predStepwise, decreasing = T),]
lm.bar <- barplot(as.matrix(bar.df), 
                  main="Important features selected by stepwise linear selection",
                  ylab = "Feature importance",
                  beside = T, ylim=c(0,12),
                  names.arg = rownames(predStepwise)[order(predStepwise, decreasing = T)], 
                  las=2,
                  col=clrs(21),
                  cex.names = 0.8)

Let us now evaluate the performance of the stepwise linear model on the test data.

lm.step.pred <- predict(lakes_step, lake_test)

# Models performance metrics
lm.step.perf <- performance(lake_test$Chlorophyl, lm.step.pred )
lm.step.perf

##                 [,1]
## model.rmse 1.3887621
## model.cor  0.8882759
## model.mae  0.9792680
## model.r2   0.7884994

perform.table <- data.frame(lm.all=ml.perf, lm.step=lm.step.perf)

Looking at the correlation we get a reasonably good result and similar to the multiple linear regression which used all the features.

We can further check the performance of the modeling using a 10 fold cross-validation with the help of cv.lm function from DAAG package as shown below.

library(DAAG)

## 
## Attaching package: 'DAAG'

## The following object is masked from 'package:psych':
## 
##     cities

lakes.cv.lm <- cv.lm(data=lake_train, lakes_step, m=10, printit=F, seed=30)

## Warning in cv.lm(data = lake_train, lakes_step, m = 10, printit = F, seed = 30): 
## 
##  As there is >1 explanatory variable, cross-validation
##  predicted values for a fold are not a linear function
##  of corresponding overall predicted values.  Lines that
##  are shown for the different folds are approximate

lm.step.cv.perf <- performance(lakes.cv.lm$Chlorophyl, lakes.cv.lm$cvpred)
lm.step.cv.perf

##                 [,1]
## model.rmse 1.3947268
## model.cor  0.8766364
## model.mae  1.0333327
## model.r2   0.7684066

perform.table$lm.step.cv <- lm.step.cv.perf

With the 10 fold CV get even a relatively better result from the stepwise linear model with \(R^2\) now 0.7684066 and the correlation between the test values and predicted close to 0.9.

Therefore, we can concluded that multiple linear regression with the 19 features selected give a reasonable or acceptable results. However, the model performance well when the chlorophyll concentration as in the middle range. Thus, other non-linear machine learning algorthims such artificial neural network (ANN), generalized additive model (GAM), etc. might need to be explored to get imporved results.

Feature selection using Boruta

library(Boruta)

## Loading required package: ranger

set.seed(123)
lakes.boru<-Boruta(Chlorophyl~., data=lake_train, doTrace=0)
print(lakes.boru)

## Boruta performed 99 iterations in 1.272344 mins.
##  33 attributes confirmed important: ABUNDANCE, AMITOTAL,
## AMMONIA_N, ANC, ATRAZINE and 28 more;
##  8 attributes confirmed unimportant: AREA_HA, LitCvrQc3OE,
## LitRipCvrQc3OE, NITRITE_N, PERIM_KM and 3 more;
##  2 tentative attributes left: ALUMINUM, LAKE_ORIGIN;

par(mai=c(1.8,0.9,0.5,0.2))
plot(lakes.boru, xlab="", xaxt="n")

lz<-lapply(1:ncol(lakes.boru$ImpHistory), function(i)
lakes.boru$ImpHistory[is.finite(lakes.boru$ImpHistory[, i]), i])

names(lz)<-colnames(lakes.boru$ImpHistory)
lb<-sort(sapply(lz, median))
axis(side=1, las=2, labels=names(lb), at=1:ncol(lakes.boru$ImpHistory), cex.axis=0.8, font = 1)
title("Feature importance as predicted by Boruta", col.main="darkgreen")

Boruta predicts 32 of the features out of 43 as important, which did not effectively reduce the features as compared to the stepwise selection. This might be due to the Boruta’s failure to narrow it down, or it could also show that problem is really high dimensional and most of the features are important.

getConfirmedFormula(lakes.boru)

## Chlorophyl ~ AMMONIA_N + ANC + CALCIUM + CHLORIDE + COLOR + COND + 
##     DOC + MAGNESIUM + NITRATE_N + NTL + PH + POTASSIUM + PTL + 
##     SILICA + SODIUM + SULFATE + TURB + OXYGEN + TEMPERATURE + 
##     DEPTH + CONDUCTIVITY + SECCHI + ELEVATION + ABUNDANCE + INDEX_LAT_DD + 
##     INDEX_LON_DD + MMI_BENT_NLA12 + MMI_ZOOP_NLA6 + AMITOTAL + 
##     MICX + ATRAZINE + METHYLMERCURY + TOTALHG
## <environment: 0x0000000015b44e20>

predBoruta <- getSelectedAttributes(lakes.boru, withTentative = F)

Let us also compare the results of Boruta with the stepwise feature selection.

intersect(predBoruta, rownames(predStepwise))

##  [1] "CALCIUM"        "CHLORIDE"       "DOC"            "NITRATE_N"     
##  [5] "NTL"            "PTL"            "SILICA"         "SODIUM"        
##  [9] "TURB"           "SECCHI"         "ELEVATION"      "ABUNDANCE"     
## [13] "MMI_BENT_NLA12" "AMITOTAL"

Here we can see that there is agreement with stepwise linear selection and boruta with 15 features commonly selected as important.

Next, we will use LASSO regularization to do the feature selection and compare it with the previous two methods.

Feature selection using LASSO

#### 10-fold cross validation ####
# LASSO
library("glmnet")

## Loading required package: foreach

## Loaded glmnet 2.0-16

library(doParallel)

## Loading required package: iterators

registerDoParallel(3)
set.seed(1234)  # set seed 
# (10-fold) cross validation for the LASSO

# Select the predictor features data only
xTrain <- data.matrix(select(lake_train,-Chlorophyl))

# Binarize the categorical features
urban <- ifelse(lake_train$URBAN == "Yes", 1,0)
lake_origin <- ifelse(lake_train$LAKE_ORIGIN == "Natural", 1,0)

# Replace the categorical data with binarized data
xTrain[,42:43] <- cbind(urban,lake_origin)
dim(xTrain)

## [1] 962  43

#Do the same on the test data
xTest <- data.matrix(select(lake_test,-Chlorophyl))
urban <- ifelse(lake_test$URBAN == "Yes", 1,0)
lake_origin <- ifelse(lake_test$LAKE_ORIGIN == "Natural", 1,0)
xTest[,42:43] <- cbind(urban,lake_origin)
dim(xTest)

## [1] 168  43

yTrain <- lake_train$Chlorophyl
yTest <- lake_test$Chlorophyl

cvLASSO <- cv.glmnet(data.matrix(xTrain), yTrain, alpha = 1, parallel=TRUE)

par(mai=c(1,1,1,0.1))
plot(cvLASSO)
mtext("CV LASSO: Number of Nonzero (Active) Coefficients", side=3, line=2.5)

# Report MSE LASSO
predLASSO <-  predict(cvLASSO, s = cvLASSO$lambda.1se, newx = xTest)
LASSO.perf <- performance(yTest, predLASSO)
LASSO.perf

##                    1
## model.rmse 1.4310867
##            0.8836414
## model.mae  1.0736569
## model.r2   0.7754114

perform.table$LASSO <- LASSO.perf

Again here, we see that using LASSO resulted in a similar performance compared to the stepwise linear regression.

library(kableExtra)
# lambda that gives best performance on misclassification error
bestlambda <- cvLASSO$lambda.min

# extract parameters at each of these lambda values
coefs.bestLambda <- as.matrix(coef(cvLASSO, s = bestlambda))

# which coefficients are non-zero
coef.nonzero <- (data.frame(variables=rownames(coefs.bestLambda)[which(coefs.bestLambda > 0)], coefs=coefs.bestLambda[ which( coefs.bestLambda > 0)]))

coef.nonzero %>% kable() %>% kable_styling(bootstrap_options = "striped", full_width = F, font_size = 12)

# Compare the features selected using LASSO and stepwise
lasso.vars <- as.character(coef.nonzero$variables[-1])
lasso.vars

##  [1] "CALCIUM"        "CHLORIDE"       "COLOR"          "NITRITE_N"     
##  [5] "NTL"            "PH"             "PTL"            "TURB"          
##  [9] "OXYGEN"         "DEPTH"          "ABUNDANCE"      "PERIM_KM"      
## [13] "INDEX_LON_DD"   "MMI_ZOOP_NLA6"  "LitRipCvrQc3OE" "RVegQc3OE"     
## [17] "MICX"           "ATRAZINE"       "METHYLMERCURY"  "TOTALHG"       
## [21] "URBAN"

intersect(rownames(predStepwise), lasso.vars)

## [1] "NTL"       "ABUNDANCE" "TURB"      "PTL"       "URBAN"     "CHLORIDE" 
## [7] "CALCIUM"   "RVegQc3OE"

# Compare the features selected using LASSO and boruta
intersect(predBoruta, lasso.vars)

##  [1] "CALCIUM"       "CHLORIDE"      "COLOR"         "NTL"          
##  [5] "PH"            "PTL"           "TURB"          "OXYGEN"       
##  [9] "DEPTH"         "ABUNDANCE"     "INDEX_LON_DD"  "MMI_ZOOP_NLA6"
## [13] "MICX"          "ATRAZINE"      "METHYLMERCURY" "TOTALHG"

# Compare the features selected by LASSO, boruta and stepwise
comm_selected_vars <- intersect(intersect(predBoruta, rownames(predStepwise)), lasso.vars)
comm_selected_vars

## [1] "CALCIUM"   "CHLORIDE"  "NTL"       "PTL"       "TURB"      "ABUNDANCE"

all_selected_vars <- union(union(predBoruta, rownames(predStepwise)), lasso.vars)
all_selected_vars

##  [1] "AMMONIA_N"          "ANC"                "CALCIUM"           
##  [4] "CHLORIDE"           "COLOR"              "COND"              
##  [7] "DOC"                "MAGNESIUM"          "NITRATE_N"         
## [10] "NTL"                "PH"                 "POTASSIUM"         
## [13] "PTL"                "SILICA"             "SODIUM"            
## [16] "SULFATE"            "TURB"               "OXYGEN"            
## [19] "TEMPERATURE"        "DEPTH"              "CONDUCTIVITY"      
## [22] "SECCHI"             "ELEVATION"          "ABUNDANCE"         
## [25] "INDEX_LAT_DD"       "INDEX_LON_DD"       "MMI_BENT_NLA12"    
## [28] "MMI_ZOOP_NLA6"      "AMITOTAL"           "MICX"              
## [31] "ATRAZINE"           "METHYLMERCURY"      "TOTALHG"           
## [34] "URBAN"              "LAKE_ORIGIN"        "RVegQc3OE"         
## [37] "LAKE_ORIGINNATURAL" "NITRITE_N"          "PERIM_KM"          
## [40] "LitRipCvrQc3OE"

Feature selection using LASSO result in 21 features being selected as important, which is the name of features as the stepwise selection above. However, only 9 of the features selected by both stepwise and LASSO regression, and 14 features were selected as important both by LASSO and boruta methods. Finally, only 6 features were selected commonly by the three methods. This indicate a strong disaggreement between the feature selection methods.

Next, let us see the prediction power of a linear regression only using the 6 commonly selected features.

select.vars <- c("NTL", "PTL", "TURB","MICX", "ABUNDANCE","PH", "Chlorophyl")
select.vars <- c(comm_selected_vars, "Chlorophyl")
small.lm <- lm(as.formula(paste("Chlorophyl ~ ", paste(comm_selected_vars, collapse = " + "))), 
               data = lake_train[, select.vars])

summary(small.lm)

## 
## Call:
## lm(formula = as.formula(paste("Chlorophyl ~ ", paste(comm_selected_vars, 
##     collapse = " + "))), data = lake_train[, select.vars])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.5734 -0.8593  0.1054  0.9602  5.2000 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.47623    0.60608  -2.436 0.015045 *  
## CALCIUM     -0.22033    0.04968  -4.435 1.03e-05 ***
## CHLORIDE     0.04758    0.03416   1.393 0.163962    
## NTL          0.90841    0.08583  10.584  < 2e-16 ***
## PTL          0.25650    0.06589   3.893 0.000106 ***
## TURB         0.83515    0.05074  16.458  < 2e-16 ***
## ABUNDANCE    0.62133    0.06741   9.217  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.563 on 955 degrees of freedom
## Multiple R-squared:  0.7112, Adjusted R-squared:  0.7094 
## F-statistic:   392 on 6 and 955 DF,  p-value: < 2.2e-16

lm.small.pred <- predict(small.lm, lake_test[, select.vars])

# Calculate the RMSE
lm.small.perf <- performance(lake_test$Chlorophyl, lm.small.pred)
lm.small.perf

##                 [,1]
## model.rmse 1.5559103
## model.cor  0.8583315
## model.mae  1.1863575
## model.r2   0.7345242

perform.table$lm.small <- lm.small.perf 

Comparing models performance so far

models.perf <- data.frame (t(perform.table))
models.perf$models <- rownames(models.perf)

perf.long <- melt(models.perf, id.vars = "models" )

theme_set(theme_bw())

ggplot(perf.long, aes(x=models, y=value)) + 
  geom_point(size=5) + 
  geom_segment(aes(x=models, xend=models, y=0, yend=value)) + 
  labs(title="Models RMSE Chart") + ylab("RMSE")+
  theme(axis.text.x = element_text(size=14, angle = 60, hjust = 1),
        axis.text.y = element_text(size=14),
        axis.title=element_text(size=18,face="bold", color = "darkred"))+
  facet_wrap(~ variable)

library(fmsb)
models.perf <- data.frame (t(perform.table))
names(models.perf) <- c("RMSE", "Corr Coef", "MAE","R2")
models.perf <- rbind(rep(2,4), rep(0,4), models.perf)
col_fill <- c(rgb(0.5,0.2,0.8,0.2), rgb(0.5,0.2,0.1,0.2), rgb(0.5,0.8,0.1,0.2),rgb(0.2,0.2,0.2,0.2))
radarchart(models.perf, pfcol=col_fill)

As we can see from the results above, all the three regressions gave very similar results on all of the performance metrics chossen on the test data. In the next sections, I will try to use artificial neural network and see if we can improve the prediction better than what we obtained using the stepwise linear regression.

Prepare the data for ANN

library(reticulate)
library(keras)
## First using only the features selected by stepwise model
# Normalize training data
step_vars <- rownames(predStepwise)
ann.vars <-c("SECCHI", "NTL", "ABUNDANCE", "TURB", "ELEVATION", "PTL",  "NITRATE_N", "DOC", "SILICA", "ALUMINUM",      "URBAN", "AMITOTAL", "COND", "LAKE_ORIGIN", "CALCIUM", "PH", "RDis_IX", "MMI_ZOOP_NLA6")
#ann.train.df <- select(lake_train, ann.vars)

train_std <- scale(select(lake_train[, ann.vars], -URBAN, -LAKE_ORIGIN)) 
dim(train_std)

## [1] 962  16

col_means_train <- attr(train_std, "scaled:center") 
col_stddevs_train <- attr(train_std, "scaled:scale")

# Replace back the factors we dont want to standardize them
train_std <- data.frame(train_std,
                        URBAN=ifelse(lake_train$URBAN == "Yes", 1,0),
                        LAKE_ORIGIN =ifelse(lake_train$LAKE_ORIGIN == "Natural", 1,0)
                        )

test_std <- scale(select(lake_test[, ann.vars], -URBAN, -LAKE_ORIGIN), 
                  center = col_means_train, scale = col_stddevs_train)
dim(train_std)

## [1] 962  18

# Use means and standard deviations from training set to normalize test set

# Replace back the factors we dont want to standardize them
test_std <- data.frame(test_std, 
                       URBAN=ifelse(lake_test$URBAN == "Yes", 1,0),
                       LAKE_ORIGIN =ifelse(lake_test$LAKE_ORIGIN == "Natural", 1,0)
                       )

yTrain <- lake_train$Chlorophyl
yTest <- lake_test$Chlorophyl

Create the model

Next, I will build the ANN model with two hidden layers and one output node for the regression. For start, each hidden layers will have 64 nodes with relu activation function. The output layer, since it is regression, I will use one node will a linear activation function.

trainX <- as.matrix(train_std)  # Features selected by stepwise only)

build_model <- function() {
  
  model <- keras_model_sequential() %>%
    layer_dense(units = 64, activation = "relu", input_shape = dim(trainX)[2]) %>%
    layer_dense(units = 64,kernel_initializer = 'normal', activation = "relu") %>%
    layer_dense(units = 1, kernel_initializer = 'normal', activation = "linear")
  
  model %>% compile(
    loss = "mse",
    optimizer = optimizer_rmsprop(lr=0.0001),
    metrics = list("mean_absolute_error")
  )
  
  model
}

model <- build_model()
model %>% summary()

## ___________________________________________________________________________
## Layer (type)                     Output Shape                  Param #     
## ===========================================================================
## dense_1 (Dense)                  (None, 64)                    1216        
## ___________________________________________________________________________
## dense_2 (Dense)                  (None, 64)                    4160        
## ___________________________________________________________________________
## dense_3 (Dense)                  (None, 1)                     65          
## ===========================================================================
## Total params: 5,441
## Trainable params: 5,441
## Non-trainable params: 0
## ___________________________________________________________________________

Train the model

Next, the model will be trained for 500 epochs.

# Display training progress by printing a single dot for each completed epoch.
print_dot_callback <- callback_lambda(
  on_epoch_end = function(epoch, logs) {
    if (epoch %% 80 == 0) cat("\n")
    cat(".")
  }
)

# The patience parameter is the amount of epochs to check for improvement.
early_stop <- callback_early_stopping(monitor = "val_loss", patience = 100)

epochs <- 500
set.seed(12)

# Fit the model and store training stats
history <- model %>% fit(
  trainX,
  yTrain,
  epochs = epochs,
  validation_split = 0.2,
  verbose = 0,
  callbacks = list(print_dot_callback)
)

## 
## ................................................................................
## ................................................................................
## ................................................................................
## ................................................................................
## ................................................................................
## ................................................................................
## ....................

plot(history, metrics = "mean_absolute_error", smooth = FALSE) +
  coord_cartesian(ylim = c(0,5))+ geom_point(aes(alpha=1))+ guides(alpha=F)

min(history$metrics$val_mean_absolute_error)

## [1] 1.041282

Check performance of the ANN model on test data

c(loss, mae) %<-% (model %>% evaluate(as.matrix(test_std), yTest, verbose = 0))

paste0("Mean absolute error on test set: ", sprintf("%.2f", mae))

## [1] "Mean absolute error on test set: 1.00"

ann_test_pred <- model %>% predict(as.matrix(test_std))

ann.perf <- performance(yTest, ann_test_pred)
ann.perf

##                 [,1]
## model.rmse 1.3320329
##            0.8977825
## model.mae  1.0034594
## model.r2   0.8054256

perform.table$ann.perf <- ann.perf

We again see that the performance of the ANN model on the test model is almost exactly the same as the multiple linear regression will all the features and the stepwise linear model. Hence, for the sake of simplicity and interpretablity in this case we have to chose the linear model with stepwise selection. However, this doesn’t mean the ANN model built here is the best we can get. There are a number of hyperparameters we can playwith to see if we can improve the models performance. We can change the number of hidden layers and nodes, we can change the activation function and the learning rate, and we can also change the model structure.

Another, machine learning technique we could try is the tree based regression which I will be training next.

Kmeans clustering

First, I will calculate TSI based on the commonly used eperical formula by Carlson, and identify the trophic state of the lakes. The result will be used to compare the performance of the clustering algorithms I will be using.

Now we have the trophic state of the lakes, we can start working on the clustering. Frist, I will use kmeans to do the clusters and compare it with t-SNE.

# Scale the data
lakes_scaled <- data.frame(scale(lakes.df[,1:42]))
lakes_scaled$URBAN <- ifelse(lakes.df$URBAN =="YES",1,0)
lakes_scaled$LAKE_ORIGIN <- ifelse(lakes.df$LAKE_ORIGIN =="NATURAL",1,0)

# Run kmeans using the four features and four clusters
set.seed(3435)
lakes_clusters4<- kmeans(lakes_scaled[,c("NTL","PTL","SECCHI","Chlorophyl")], 4, nstart = 20, iter.max = 30)

Clustering using t-SNE

Next, I will do the clustering using t-SNE and compare both kmeans and t-SNE with the emperically calculated TSI results.

library(Rtsne)
set.seed(3434)
tsne_lakes <- Rtsne(lakes.df[,c("NTL","PTL","SECCHI","Chlorophyl")],
                    dims = 2, perplexity=80, verbose=TRUE, max_iter = 500)

## Performing PCA
## Read the 1130 x 4 data matrix successfully!
## OpenMP is working. 1 threads.
## Using no_dims = 2, perplexity = 80.000000, and theta = 0.500000
## Computing input similarities...
## Building tree...
## Done in 0.47 seconds (sparsity = 0.266508)!
## Learning embedding...
## Iteration 50: error is 54.431234 (50 iterations in 0.20 seconds)
## Iteration 100: error is 52.519277 (50 iterations in 0.18 seconds)
## Iteration 150: error is 52.512390 (50 iterations in 0.16 seconds)
## Iteration 200: error is 52.512236 (50 iterations in 0.16 seconds)
## Iteration 250: error is 52.514158 (50 iterations in 0.16 seconds)
## Iteration 300: error is 0.680049 (50 iterations in 0.19 seconds)
## Iteration 350: error is 0.584024 (50 iterations in 0.19 seconds)
## Iteration 400: error is 0.562367 (50 iterations in 0.20 seconds)
## Iteration 450: error is 0.555057 (50 iterations in 0.20 seconds)
## Iteration 500: error is 0.551683 (50 iterations in 0.20 seconds)
## Fitting performed in 1.83 seconds.

Let us plot the results of kmeans and t-SNE and overally them with the calculate TSI values.

clrs2 <- c("darkblue", "deepskyblue3","chartreuse4","black")
trophic.col <- trophic
levels(trophic.col) <- clrs2
trophic.col <- as.character(trophic.col)

par(mfrow=c(1,2), omi=c(.1,.3,.4,.1))
plot(tsne_lakes$Y,  t="n", sub = "t-SNE vs calculated TSI", col.sub="brown" ) 
text(tsne_lakes$Y, labels = trophic.num, col=trophic.col) #, pch = lakes_clusters4$cluster)
legend("topleft", levels(trophic), text.col = clrs2, bg='gray90', cex=0.8, pch=levels(trophic.num))

plot(tsne_lakes$Y, pch=levels(factor(lakes_clusters4$cluster)), 
     col=clrs2[trophic.num], sub = "kmeans vs calculated TSI", col.sub="brown" )
legend("topleft", levels(trophic), text.col = clrs2, bg='gray90', cex=0.8, pch=levels(trophic.num))

mtext("t-SNE vs K-means Clusters for US lakes", side = 3, cex=1.5, outer = T)

# Check k-means accuracy
caret::confusionMatrix(factor(lakes_clusters4$cluster),trophic.num)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   1   2   3   4
##          1   2 173 214   0
##          2   0   0 241 119
##          3 151  52   0   0
##          4   0   0   1 177
## 
## Overall Statistics
##                                          
##                Accuracy : 0.1584         
##                  95% CI : (0.1376, 0.181)
##     No Information Rate : 0.4035         
##     P-Value [Acc > NIR] : 1              
##                                          
##                   Kappa : -0.0843        
##  Mcnemar's Test P-Value : NA             
## 
## Statistics by Class:
## 
##                      Class: 1 Class: 2 Class: 3 Class: 4
## Sensitivity          0.013072   0.0000   0.0000   0.5980
## Specificity          0.603889   0.6022   0.6988   0.9988
## Pos Pred Value       0.005141   0.0000   0.0000   0.9944
## Neg Pred Value       0.796221   0.7078   0.5081   0.8750
## Prevalence           0.135398   0.1991   0.4035   0.2619
## Detection Rate       0.001770   0.0000   0.0000   0.1566
## Detection Prevalence 0.344248   0.3186   0.1796   0.1575
## Balanced Accuracy    0.308481   0.3011   0.3494   0.7984

We see here that on the left figure, t-SNE has a strong agreement with the calculated TSI. The figure on the right is a bit complex, the location of the observations are based on t-SNE coordinate, the colors are based on the four trophic stated as calculated by the TSI, and the labels(numbers) are based on kmeans clustering. Thus, we observe that kmeans (also from the confusion matrix) has a very poor agreement with both TSI and t-SNE results.

set.seed(33)
tsne_lakes3d <- Rtsne(lakes.df[,c("NTL","PTL","SECCHI","Chlorophyl")],
                    dims = 3, perplexity=80, verbose=TRUE, max_iter = 500)

## Performing PCA
## Read the 1130 x 4 data matrix successfully!
## OpenMP is working. 1 threads.
## Using no_dims = 3, perplexity = 80.000000, and theta = 0.500000
## Computing input similarities...
## Building tree...
## Done in 0.49 seconds (sparsity = 0.266508)!
## Learning embedding...
## Iteration 50: error is 54.563281 (50 iterations in 0.45 seconds)
## Iteration 100: error is 52.482003 (50 iterations in 0.28 seconds)
## Iteration 150: error is 52.475115 (50 iterations in 0.25 seconds)
## Iteration 200: error is 52.474545 (50 iterations in 0.25 seconds)
## Iteration 250: error is 52.473286 (50 iterations in 0.25 seconds)
## Iteration 300: error is 0.621888 (50 iterations in 0.34 seconds)
## Iteration 350: error is 0.509386 (50 iterations in 0.38 seconds)
## Iteration 400: error is 0.483002 (50 iterations in 0.38 seconds)
## Iteration 450: error is 0.473704 (50 iterations in 0.40 seconds)
## Iteration 500: error is 0.466167 (50 iterations in 0.40 seconds)
## Fitting performed in 3.38 seconds.

tsne3d <- data.frame(tsne_lakes3d$Y, colr = trophic.col)

info <- left_join(select_lakes2[, c("UID","URBAN","LAKE_ORIGIN")], lakes_info12[,c("UID", "NARS_NAME","STATE")])

plot_ly(tsne3d, x=~X1, y=~X2, z=~X3, colors = trophic.col) %>%
  add_markers(text= ~paste('Lake name:', info$NARS_NAME, '<br> State:', info$STATE)) %>%
  layout(scene = list(xaxis = list(title = 'D1'),
                     yaxis = list(title = 'D2'),
                     zaxis = list(title = 'D3')))

As it can be seen from the plot above t-SNE was able to cluster the lakes in the same way as the emperically calcualted TSI levels, however kmeans had a very low performance with prediction accuracy less the 50%.

Discover clusters using all features

In this section, I will try to use the entire features and explore is there is another clustering among lakes. I will start with kmeans first, and then use t-SNE to answer this question.

library(cluster)
k=seq(2,100,4)
lakes_clusters=vector("list", length(k))
within.ss=vector("list", length(k))
between.ss=vector("list",length(k))
sil=vector("list",length(k))

set.seed(3456)
for (i in 1:length(k)) {
  lakes_clusters[[i]]<- kmeans(lakes_scaled, k[i], nstart = 20, iter.max = 40)
  within.ss[[i]] <- lakes_clusters[[i]]$tot.withinss
  between.ss[[i]] <- lakes_clusters[[i]]$betweenss
  
}

plot(k, within.ss, xlab = "k values", ylim = c(1000,50000), ylab = "sum of squared errors", type="b", lwd=2, pch=18)
lines(k, between.ss, col="royalblue2", type="b", lwd=2, pch=20)
legend("bottomright",c("Total within clusters","Between clusters"), lwd=c(2,2), col=c("black","royalblue2"), pch=c(18,20))
title("Clusters Variance vs Cluster Numbers")

clusplot(lakes.df, lakes_clusters[[1]]$cluster, color=TRUE, 
         shade=TRUE, labels=2, lines=0,
         main = "kmeans with k=2 clusters")

set.seed(11)
tsne_lakes2 <- Rtsne(lakes.df, dims = 2, perplexity=80, verbose=TRUE, max_iter = 500)

## Performing PCA
## Read the 1130 x 45 data matrix successfully!
## OpenMP is working. 1 threads.
## Using no_dims = 2, perplexity = 80.000000, and theta = 0.500000
## Computing input similarities...
## Building tree...
## Done in 0.51 seconds (sparsity = 0.267741)!
## Learning embedding...
## Iteration 50: error is 60.389727 (50 iterations in 0.21 seconds)
## Iteration 100: error is 56.506053 (50 iterations in 0.19 seconds)
## Iteration 150: error is 56.492306 (50 iterations in 0.20 seconds)
## Iteration 200: error is 56.463207 (50 iterations in 0.20 seconds)
## Iteration 250: error is 56.452528 (50 iterations in 0.20 seconds)
## Iteration 300: error is 0.820401 (50 iterations in 0.17 seconds)
## Iteration 350: error is 0.746273 (50 iterations in 0.17 seconds)
## Iteration 400: error is 0.727430 (50 iterations in 0.17 seconds)
## Iteration 450: error is 0.719992 (50 iterations in 0.17 seconds)
## Iteration 500: error is 0.715898 (50 iterations in 0.17 seconds)
## Fitting performed in 1.86 seconds.

plot(tsne_lakes2$Y, main="t-SNE Clusters for US lakes", col=lakes.df$URBAN)
legend("topright",c("Urban", "Rular"), col=c("black","red"), pch = c(1,1))

Hierarchical Clustering

hc.complete = hclust (dist(lakes_scaled), method ="complete")

plot(hc.complete ,main =" Complete Linkage ", xlab ="", sub ="", cex =.6, labels=info$STATE)
rect.hclust(hc.complete, k = 4, border = 2:5)

# Cut tree into 4 groups
sub_grp <- cutree(hc.complete, k = 4)
table(sub_grp)

## sub_grp
##   1   2   3   4 
## 559 517  53   1

We notice here that cluster 4 contains only 1 lake and cluster 3 contains 53. The majority of the lakes are found in the first two clusters. From this we could say we have two major clusters.

We can also use the function fviz_cluster from the factoextra package to get a 2d visualization of the four clusters.

library(factoextra)
fviz_cluster(list(data = lakes_scaled[,-43], cluster = sub_grp)) # URBAN had 0 variance, remove it

Cluster 1 and 3 seem more distinct when cluster 2 could be split into these two clusters. One observation that falls into cluster 4 could anommaly or could be in cluster 3 if we see it from low dimension only.

fviz_nbclust(lakes_scaled, FUN = hcut, method = "wss")

fviz_nbclust(lakes_scaled, FUN = hcut, method = "silhouette")

sub_2grp <- cutree(hc.complete, k = 2)
fviz_cluster(list(data = lakes_scaled[,-43], cluster = sub_2grp))

table(sub_2grp)

## sub_2grp
##    1    2 
## 1076   54

hc.wd = hclust (dist(lakes_scaled), method ="ward.D")
plot(hc.wd ,main ="ward.D method ", xlab ="", sub ="", cex =.6)
rect.hclust(hc.wd, k = 2, border = 2:5)

sub_2g.wd <- cutree(hc.wd, 2)
table(sub_2g.wd)

## sub_2g.wd
##   1   2 
## 742 388

fviz_cluster(list(data = lakes_scaled[,1:42], cluster = sub_2g.wd))

par(mfrow=c(1,2), oma=c(1,1,2,0.2))
plot(tsne_lakes2$Y, col=sub_2g.wd, sub="With ward.D method")
legend("topright",c("Cluster1", "Cluster2"), col=c("black","red"), pch = c(1,1), title = "Hierarchical clustering")

plot(tsne_lakes2$Y, col=sub_2grp, sub="With complete linkage method")
legend("topright",c("Cluster1", "Cluster2"), col=c("black","red"), pch = c(1,1), title = "Hierarchical clustering")

mtext("t-SNE vs Hierarchical Clustering of US Lakes for the Year 2012",side = 3, outer = T, cex = 1.5)